{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## EM algorithm to solve Mixture of Gaussian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-24T04:15:15.560386Z",
     "start_time": "2020-03-24T04:15:15.556146Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-24T04:34:47.290481Z",
     "start_time": "2020-03-24T04:34:47.122466Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.01501703 -1.96645385]\n",
      " [ 5.26279743  2.78817091]\n",
      " [ 7.80861827  5.11821458]\n",
      " ...\n",
      " [-6.28980122 -5.54954711]\n",
      " [-0.6126305   0.02192317]\n",
      " [-8.2015517  -5.49554098]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1a26326bd0>"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAD5CAYAAADGMZVsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztvXuQHMd54PnLKgw5g9dgQEwNAJJ4UEMIINUUogk+JFLL9VoyBenW8tpaXsi0T3LYOq/D67UvrGBsyOFQnC8oO+QNrfZky/bKsbZ2RduitSdbIQvUUjZNk4RkiWxSbIkYAqCIl2Yw3eA88JgZAlOV90d2ddcjq7r6NdM9nb8IxKDrkZWVVfXll19++X1CSonBYDAY1jbWalfAYDAYDJ3HCHuDwWDoA4ywNxgMhj7ACHuDwWDoA4ywNxgMhj5g3WpXQMenPv8NAdwEXFztuhgMBkOPsRk498hHHwy5WrYk7AvF0sPAn0Q2bwB+K59zPhk59vvAXsCvwFP5nPOvE4q+CTjTSt0MBoOhj9kFnA1uaEnY53POY8Bj/u9CsfR/AL8F/GHwuEKxNAjsA7blc04Wbf0iwJG/+SLLy9daqWIMIQQ7d48zefokZo2BwrRJGNMeYUx7hOnm9li3boDDP/VzoLGKtM2MUyiWdgD/L/DefM6Zj+x+G3Auo6Cvsrx8jeVrV9tVRQCEsHBdl+Vr15DSa2vZvYppkzCmPcKY9gjTq+3RTpv9/wP8z3zO+bZm30FAFoql7wB7gGeBf5/POZNpBQohEKK9c8h+ee0ut5cxbRLGtEcY0x5hurk9hBCJ+9oi7AvF0nbgQygNPonvAh8D5oHPAF8C3pVW7s7d47iu244qxsveM96RcnsZ0yZhTHuEMe0Rphvbw7btxH3t0uw/BPxjPue8rtuZzzl/Cvyp/7tQLD0CXCgUS1vzOWcmqdDJ0ydZvtZum73Fzj3jTJ462VNDsE5i2iSMaY8wpj3CdHN7rBsYIH/ofv2+Nl3jXwN/nrSzUCz9IvB6Puf8Q2XT9SivnDfTCpVSdqwxpfS67kGtNqZNwpj2CGPaI0w3tkfahHHLwr5QLFnAXcAvpxw2BvyHQrF0GLgEfBr4Sj7nXGn1+gaDofvxhsZgyIHFEtbi9GpXpy9ph2Z/A7ARmApuLBRLR4BnKv72nwIc4EVgEDgC/FIbrm0wGLocd+cDMHonWDZ4Lm75BezJp1e7Wn1Hy8I+n3PKQGwKOJ9zDgf+vwz8RuWfwWBY4/iavPTcmqAH9Xf0TrzZCazFaaPxryBdGS7BYDD0LlFNvirofSwbhhzckf1G419Bus9R1GAw9Cze0Fhck49OGnou0lvWa/xDYytb4T7CCHuDwdA+hpy4Jh9c6CM9KD2PsNYlavyGzmDMOAaDoX0slvSmGx8pEXOvqv9Hj/Ncdb6hIxjN3mAwtA1rcRrKLygNXnuA0t6rx3mVFfKeC6XnzSRtBzGavcFgaCtidgI5eqfGR4+Q9m5PPo03O5HojWM8ddqLEfYGg6G96Oz2oNXercVp0Ahy45vffowZx2DoMN7QGN7WXP94mvh2+yDSg9NHEHOv1m0LrUeP8dRpGaPZGwwdpB81VGtxGrf8QtjXvvQ8rHeQuw/XbwvdyMD31GnAnGPMQGGMsDcYOkSihlpZPbqWidrjAeS+h7O1hc6jp0FPnX7sZOthzDgGQ6dI01D7AGtxGmumqIR5A23RqqeOMQPpMZq9wdAp2qChrhkabIt6njqptMkMtNYwmr3B0CF62Ze83ZPKzbRFaGTQCNoJYonceHODtV5bGM3eYOggLWmoq0Tc3v1doNxyuc22RaMTrdbiNO7sMdh6ey1UgxAwsh+v/EJPPINOYIS9wdBhknzJuxGtvXvbQbj8ZFvKb7Qtmp1oFZfPIW+IpMTWmHL6yWPHCHuDwVAjyd7dJmpx7pdVMLQUIZvVm0krsDPMEcQ6kktnEFPPrFmhb4S9wWCokSQk20BIuEqJFCKmrQcFd5aJ1iTNP8nX3xfk2o5keC9y8y7c0tp00zTC3mAwVNEKyfKLMNRaue6WA2Hh6tvSLRucQ3izE8hIMhNmj6Vq5/U0f3vyadyFEmzaA5dOYc8dq5WTFNJBrN21EEbYGwyGEGJ2QqUTlBJx8ST2Uhn27mu6PHfnA+AcApHg/Ccs5K4HYXBbWHCPHIDZCRjZr9XO62n+Ia1/6wHc9U5NY08LxbxG3TSNsDcYDFWiZhFpr4Op5j1xqtp3kqD3CQp6H8tGXD6rXDYjNnlvaAw5sBmkq7Tx6gWV5l9P64+NYEKVXptrIYywNxgMQIpZZG5Cf2wWL5Ykc0kUy1bB0oKdgufVyg9cI9wheYBUZiEpYUEd723N1bX3+66gcse7YPMu1WkERg9rzVPHCHuDwaBIMosMjgLL1U0NuUPqzCVShlMVghLa0fj3mnj48Q4p0DkIAesddUzGFbvW4jT88Msxwb4WY+u0LOwLxdLHgE8CVwOb9+VzzmTgGAv4NPDzgAd8Jp9zHm312gaDoY0kCcilMjCifjYY3C1mLolq7z5CxDsAYcVt5/VGCn4mrJliqjeOrp6keeqsgUnbdmj2B4HfzOecz6Yc82vAvcA+1FvzZKFY+k4+57RnpYbBYGiJqmarnRCtCfsk7V9uHk+c0Kx6xYzdA0Pb9BUQIm5/ly5yYBPe0FjYfz56XOhGsmfCSmSNxtZpl7D/r3WO+Vng0/mc8wbwRqFY+hzwS0CqsBdCIOpN7DSIX167y+1lTJuE6bf2cHfcp1bJWjbCc+Hia7C8BJdPY88dD7WHWCqBuxQXhtvvxLMt7KnntOWLSvmwDDKhIp4HYllp9NJTx+24G8buxLvwEvbUc9hLZdyLr8Pm3fHzpYSZY9hLZbz1YzA4ilgqY83+QO3P+Dy19+i5IK+pNuji90NER0fBfVImtXx9CsXSIHAJ+AZwNzAFfDyfc/4uctw88C/yOed7ld8/BfxOPufcoSv3U5//xjAwV3j+WVy3PQs6DAaDYa1j2zb5Q/cDbHnkow/OB/e1qtk7wFHgs8BTwE8AXyoUS3flc05gBQMbgMXA7wVgfb3CJ0+fZPnatRarGEYIi517xpk8dRIpvbaW3auYNgnTT+3hjdyG3PXulANcROm73Lx+LtQe7viHYMNo/HjfJi9dKL+EWJpNL78BxJlvYs2+oq6/9wN67X7pEly/MWz/9yeEPRcqI4QsuFveCrveE/EQcrFee5ybxka68v1YNzDgC/v4vlYKzuecM8ADgU1fKxRLTwGHgaCwXyC8Bm89cLle+VLKjjWmlF7XPajVxrRJmH5oD7lYRkoredLTHoDt98LFJ1je/g6syWfUatj1Y3q7uS9jxQA4dyNPfR3Sys9cUYncuAsx8331e+oocvMtcdPM4IhmojdwL6N34c4ez2a/F+uQ1vXhbfYA8nrlndSN70eapaYlo1OhWHp7oVh6JLL5emApsm0CNTnrs6+yzWAwrCKxOPNpZt1tB3F3HYY978smvIUFgzfE49hfqUyyNoIforgSX99anIaLp/XHpdFIpjBdXPyqd1Lv0aoZ5yLwiUKxNAF8DfgZlNfNhyPH/SXwSKFYehrYBPwK8OstXttgMLSBoNeKdO5Uq1l1QtMPYaCbmNT5zgNIDzF7AukugxCI+RPKHfOWD8Lw3sYqatnI4VurHjFi6hnkpl31ffgj9ZQbb4KZYv3LJQZTC3gn9RCtmnFeLxRLHwJ+FyXQTwA/mc85U4Vi6QfAJ/M55zGUTf9G4GXUoOrT+Zzz1daqbjAYWiG6kMgDuH5rsrCUHlgDmoJcuDQJm24ML3KSHmzeg9x+by38gmXD4rReUGdh7G5cy65Ftpw9pjogXxgvTMOGBBMTVEYIB/DKhUymHK37Zhd64WShZdfLitCOCe58zrk98P9l4GOVfwaDYZXRrRAVSzPpwnf2OIzcFtGkPbj4OmzeqwR9SLMWsPHG2rGRxUnhxVZJI4PI9kAZcmR/TdBLDxbOI879PXI0D1tvSxbKDfrMr0TymZUIzdCbXZTBYGiapBWict2GuM1eSigrk4d99n/F7e9vvKLMMdHQxdH/+1g28sYfDxxTvVAlzk3k2u5VYlg2cng8EjLZgo03Ivc9HBf00XvyXOR1w23Lr9sq7s4HkPseRu5+L3Lfw6oj7gAmNo7B0G8krRC9XufJImB5Aa4Hb2g0ZtaQO96VbDJJYuNO3L0/DZv31M4VFrHVVkKAfV38fM8FRHJ44ihC1FxCPU/93vEO5NjdsZg3Kx38bCVDMxjN3mDoN5K8TC6d0mz3YPshAOT4Q7g7H1Ahgv0Jzk27Gr++EOHRQHC77tgoF08j5k/ERwJJeC6c+jpMHQUha1q/L1grGv5Kadgh0kIztBkj7A2GPiPmblnxMrHnjsW3C2rat2WDcyfu2DtrGrBOk/aFcJobpx8SoR4xs5KHOP9spYwM5wXuTVy9GB+FVARroobdaVNPUsfbgXj6xoxjMPQhSUHCQm6Y1w3DjneEfeKFDTvvUx42urSBoITw0gxcPwwkmHikhPkfwrBmYVQavuwfctLPEwIufB9xoRAOopYQ9lgO31o3+JnfwYml9gnierly24kR9gZDn1J1txxy8Cq//e1UknfIsbv1GrQubaCPsGBwa/rFhYCrV2DuNdjylpSUhZGLWxZyW16ZnNKiX0ovLOgr9+UulGDD9lCyk2ru21gZEumpOP4h7yV3CS63L2Bv09E5G8QIe4Ohh2jnBGK9BB1VrXObNl6hOu/qReV6uWW88QqM3lF/xWsUKWHb25THTXR7KB5O/FRvaAzWO7XjKslO2OAkhH4QCGud3sSDmrAWC+0RzCvh3mmEvcHQI7Qze1JWLxB78mnk3ARsH1EmBjuwqEp6MHZ3up98kOgxWQV98Dz/rxUZCehGADveBT/8cm1b0mRoEr7tPGluYnBULeLqEcwErcHQA7RrAtEbGsPbmku3UScxdzziARNwf0xcdRtQsRvV4v3zmzkPYPPucPskTYbq4vTIgO1cdx70XIwcI+wNhl6gDS56QddCxu6Oe8NovEDcnQ8gxx9SP7bsC9vvswjhZgW1z+zJdK+e1GtbofZJ8kJSgdkq15ASluYQrz6GPfVPyedBJUZO72DMOAZDL5AxgXYS2pGB5wGVSU6NF0j1HOHVzknzgNFp4TFbegOauuciSt9CXrexNqna4LX9CVaf6GQogHTuCpuJrtuMvH5ryIZeTa24aQ9cfr0X46AZYW8w9AItu+hpRwYWTB5FXLuon/D1z8niD++5MHMMBjaohVbBwGTrx2q/swY+k17t/k58UcXQv+HtsOmm7IuvhIAbH4C5Y6HNwclQd/t9cfu/ZcGe9+GWnOqcSGi+ZGS8rd44K4UR9gZDj9CSi17CyEBcPJlcjn9OPUVcejB1FLv0beWuuXCr0rTfnFFRNBcvwPIVWLdBeeBkIqyp23PH8Kx1yM03Zzy/wsBG3L0/jTj/XEK7JdycsKoT1kDHvXFWAiPsDYYeolkXvaSRAYC7/X6QMib467peQsV8YsGOd+IOjytXRn80IKlpzdKDmVdqMWrqITTeQYslZXqKauKp5ajQDHJ4T9VcFfRiEvMnkNvv0dcpOCeyBrxxjLA3GPqEWBCzkf3It/5cVdDJ7ffilp8PuXOmul5CwBXSDtvVhRWZzLVg64HGVs3qQhE3M98bvFbExdRanMYtPQ/OnXFf++CciM4EZbxxDAZDtxIKYubcGRGEltads+p1cuElvQuiT72JV2E3Fh7Bc8MTrGkhEoLeNPWoZLzytubwhsawJ59GvPoYzL9ec8MMzIkYbxyDwdC7DCWsGk1J7GFPPceytT7Z7p4hJWBmQS+lqsuew7gXb4fL55QnTNI1oouu0upVWQwmo4vTfvjlxBXKwVGRWCqpkU6PYYS9wdCPLJb0sWVS3DndHffBDbdr91Vt9GkpDRsJeFYV3rYKh7x5T/M++0LUzDDSBSnCk60Rs07SnIi/TwiLXvS9NGYcg6EPsRanofRC2K3S89LdObcdjNutpQeXf1TRxBMSkV/+kYonr1upmpVWFmd5Lpw+gjj9BEz9s8bVsrHFad7QaOhvr2A0e4OhT/FNE3L4Vq03TgydR0q5qLYHc836zJ1EnD+qJkF3PqA0al9mJy22ki5gZRPuQY+fxBy2gVj9qIVicvs9TS9Oc3c+gNh2B1x+Ejn+EN6Fl5uOT7TSGGFvMPQxDblyShn2hpGyFrkyKmw9tyroa6t3gxq1rJh2KqaVN44hrpxDDmyGne+sX4/5M7A4CUtvwOANMHI7DG6OHzv1z9jTz4Xut9nFadoVxR1KIdgJWhb2hWLpA8AngZuB48Bv5HPOs5FjtgElYCGw+XfyOedTrV7fYDC0j9QQyrr8tMH/+wI/KkB1q3eDmaok4C5gzRT1mncUIWD4ZtiyO127rywaixJK0OItV8MY1xXYuhXFKRPa3UZLwr5QLN0C/HfgJ4FngA8BXy0US3vzOWc+cOhB4MV8ztFkCDAYDN1A3RDKOj/7IELAQhlx5khIcErPTTDbRHLBVjTkkOatI3hu0J8/rcOJYC1O41aSlsS8cgIEOz/tiuIOpRDsBK1q9ruBz+dzjt9CjxWKpf8C7Af+OXDcQeB7LV7LYDB0iLT49ra/eOhKCTbXSTA+uBV5/Va8wOhAWDYybVTgX6+iIVfnEjaPw457wh5D9dw7dekIG7xf/7xo50f5BfXPX1HsuVBOv0430ZKwz+ecp4Cn/N+FYuleYCPKnBPkIPCWQrF0AhgC/gr4eD7nXE0rXwhRcXNqH3557S63lzFtEqYf20MMjSpbdNBEIdR28eYbAFjn/xFv4/9eJ/crsOvdVQHpXXgJMXdcpfJLM81ID7FUqra5vVSGpTKubdW8gKQHl8/DhrHksjwX8UYBa6mcWs/U+10qq7g32+6oHSOAbXcgTj6OmJ+AsRGs1x5HLKRfZ6URKR2hkM3Gio5QKJbGgaeBz+Rzzu9H9v0hMAv8PrAF+DLwRD7n/LaurE99/hvDwFzh+Wdx3RbctQwGg6GPsG2b/KH7AbY88tEHg6b09njjFIqle4CvAp+LCnqAfM751cDP+UKx9HvAJwCtsPeZPH2S5WvX2lHFKkJY7NwzzuSpk8gsoVv7ANMmYfq1Pdwd99W0aM+F8ovY54/G2sMbGkXecBBuOJCpXHHmm8jBERg9WPG+CUyqei7MHsc+pw8Z7A2NquQpUVfJS6crC60s5c1z6VzF+yd7CIPY/c6dQFw+V415o7uuOPk49tIbXft+rBsY8IV9fF+rhReKpfcBfwn8Zj7n/KlmvwAeBf4on3POVjZfDyzVK1tK2bHGlNLruge12pg2CdNv7WFNPoM3ezzkjRMc9/vtIRamkfIl5Nbb68en91yk58HoXYEUhigBPfVtxMXXYtcJIgcdpD0Y3mgPwPCtgUnaAdi0F6aONvS8gvcrN94EI7chb8jV7PMXXo7Y7AtYC9NIP3BcF74faZaaVr1xdgNfAj6Szzn/U3dMPufIQrF0CHi0UCz9MjAGfBz4k1aubTD0OqlujqtEVr/7bF4zakWusGzl8RJE2Ihrl+rfd1KGrqQUjQ22o7U4jQdw83tik7Xi+GPQbP6ALqRVzf7/AjYAXygUS18IbP/fgK8Bh/M55xngF4DPAZPAVZSg/4MWr20w9Cx13Rx7gNAK3LG74wL5jR8g5l6t/W5i1aq1OI07ewxGDtQ07JljKlxypDx53XA2f/koKfl9rZliT/jQZ6FVb5zfAH4jYffGwHE/Aj7QyrUMhrVCFre/XsEfCbiVe6jmthUCRu9A3nB7zWWxzqpV3UjH3flATdBLD2YnsM8ewXUX4tfb8Q7k2N2Nd5wt5vftFUy4BINhpUnRJDupRXbSbFTV8m/6V7DhxnBSkwwmEd1IR8xOhDtFYcHIfryKMFejinHYHvDFb6LjbDm/b49ghL3BsNKsgibp7jocMoV0zGy0YWfigqkkk0jSSEfWsc2ruDsOMhqmuYmOs6X8vj1C96wGMBj6BG3mow5qku6u90HQc8bXfiMZqXy8obFqFqeGSMokVenIEstNGulIGc+MFe0U/Y4z7ZiM+Fm81qKgB6PZGwyrwkppkt7QmJrMTAlPEERnTlk39Uy2i+lGLFLCzCvIShwa7cgiYaQjLp5E2utSzSv9YoJpB0bYGwyrREPhhZvAGxpDbsunatvR47UTx3MTma4XE7zSg5ljiAsvIvc9nJodKlFgL07X7RT7wQTTDoywNxjWICENXRdxcuaVuFBMMqcMjgLLZEEneL2tuboT0mkCO0un2OmOcy1ghL3BsMaIaejB0L8Vbds++0T8xKSJ46UyjeRcjQnejBPSRmB3FjNBazCsNbTJQiqhf1/9IvaZr2tPS544zh5vprFyjWBfSYxmbzCsNZImPOvEeIcEc0obQvgau/rqY4S9wbDGaNVDpVPmFGOmWV2MsDcY1hD+KlkxO7GmgngZWscIe4NhjaBLo9drwdUMncNM0BoMa4BEH/lGV8Fqyx4N/V1Nml7dazCavcEQxBsaU/lJE/zKOxVMrOVyOxRczd35gMrFevlJ5PhDeBdeVpOtqxCLP2tY6G7ME9ANGGFvMFSoChPhwcUncHfchzX5THx/m4OJZSm3rgDrQHC16mhBVLIxVUYL7rr1KxNUTVeXOmGhE6NnGuFvhL3BABFh4qea23YQb/Z4JbpiZ2LQZyk3S2fQiAdOZs3XHy0EU+9Zdk3Qt7Ed6tYrw8hF35aHkM6dKgRyjyaJaRdG2BsMUF+YNGkmqStY65TbSCeTxZe9odGJP1oIxlCTXsfMRan1yjJy0bZlYFqyh5PEtAMzQWswQP1QuU2E0nV3PoDc9zBy93uR+x5WAq3R66Z1BhrSwvQ2OomrXfk6c0xbX+ktNz1xmqVemVbh6toydlPJbbfWMZq9wUDEDOJrsuUXq8Kk0YVKWTVya3Ead6EEG7bXYtgsBLTydtrimxid2JNPI+cmYPsI4uTjWAvTuMtXwu2wMA27D6uk4s2YSjLWq97IRfuMhAivAF6D6QazYoS9wVDBFyZiaBS2gH3+KFKzvyF7dxCNAPOGxmC9U4tKKQSsd6qJs9sar73JjkPFxhnBWiwjUe3gLpRg0x54cw52vKM1G34D9aq3Cjf6jKJx9Ps5Jo8R9oauoFvc5azFacRSGbbsS9yfyTadVYBl6BTaFVemXR1HOHyyF4+d06ANv1avQzUbuxBKUDdxr6FnlCEefr9ghL1h1emUS+Nqklmwdjj8b7QTbbXjiIdPtuLx8pswlYjZCeU1U91gtW0y1cTkUbQs7AvF0j3AHwP7gO8BH87nnBORYyzg08DPAx7wmXzOebTVaxt6n065NKZeb4W0vCyCtV0uk7p9SZ1oS8IvKXyy32E1ayoZcpR7ZJA2ePkYarQk7AvF0iDwFeBjwJeB/wj8OXBf5NBfA+5FdQgjwJOFYuk7+ZzzZCvXN6wBOrTyU8dqjCCyCNZWXSYTFxI12YmmdSrSW46bbjwXTh9BWOua6kTdLQdg87i+3D6dTO0ErWr2PwbM53POXwAUiqVHgd8sFEsH8jnnWOC4nwU+nc85bwBvFIqlzwG/BBhh3+90YOWnjpUeQTRKWqeQVndAu09G29TfV6cTTUs47u64D0bvCptuKpq8PXcsscw03Ft/LuyJFCm3G57NWqFVYb8fqGYjzucct1AsvQ4cAI4lHQecQJl0UhFCINqQOCFcphX6a1jdNrGXyrjl78K2g4FojS9iL5XbkjTDRwyNqmX/wdWgQm0XS+FMTPXawxsaVXlZl8otZ3HKUl5a3dX/4/uQLrhLsU5ULJVS70tsu6NWngC23YGcV5+uuOH28LVcF84+iT13vKln5W55K2L9DcAyVbcnKeHCK4iZl1VbdOF32s0yRATnTiK0Kuw3AIuRbQvA+jrH6Y6JsXP3OK5bZ5FEk+zcM96RcnuZ1WuTMlwODPKGgL16b5jmWYaLmryrW0j0vElvj2WURTJ7btZ00spLqTvo920g3KY+2+vUWXfOmDr+5iv/GN83Aow0+6wkXPpGfPN11K9nF9CNMsS27cR9rQr7BdSnGWQ9cLnOcbpjYkyePsnytWstVTCKEBY794wzeeokMqgN9TH90ibujvviI4jzR2PHJbWHNzSKHH8ori2ffLwpDT+pPOaOqw7Ir+eFl9S+hLqn3Vcjo5DE+zv7JDePwNmN/xJpDbbl3gHcLftg14NhTx4pAak0+sq921PPNVV+p+jm72XdwAD5Q/fr97VY9gTwi/6PQrFkA7cQNtn4x/neOlT+Hz0mhpSyY40ppdd1D2q1WettYk0+gzd7PDT5KFOOj7aHHHSQ9mD4IHsABh3kQuO25cTyRm6vCVx7AEbvQhx/DBLqnnZfYmFarXCF1Hv1j/UWZ8M2dO8q3PweFeKY65DYNUFcLmAtpLdhGtbsBO62Q+HrCWrCv3LvbiUYXbfRjd+LlMlPo1Vh/xRwQ6FY+gjwFyhvnBP5nBMV5H8JPFIolp4GNgG/Avx6i9c2GBqmJbfDdk8m68pLCDQmN49jTz+XWPd2+JJrV/OuW081tr9VmZg9/y3E/Im2CGD7xBeVN86mPereR+8IH2DcL9tGSzMM+ZyzCLwf+FXgDeA9wEMAhWLpB4Vi6eHKoZ8FngZeBo4Cf5zPOV9t5doGQzvJkgEpUzCuBsgcaAxgxz36QGoN4A2N4W6/H3fsPv19JvnQhyptI67Oa/38mw2EZs8dwz57BPHGSw0HmzNkp+VFVfmc8wJwl2b77YH/L6N88T/W6vUMhnbj7jocScbxXUBvh25kBWqWBVxidgLpLiuhuvQGwlqHRMLIbeHwvCLdVbTetdydD4BzqOrdIrffi1t+vn4Y4diFahEu/Wu1a/1CW+MAGWKYcAmGvsbd9T7YeltNg7VsNdmp80qpkMVkkkUAxmLMSJCWpdwmtRfWmzTqXcsbGgPnzrAboxUPRxATttEwCFLC0kw4wuXsBIzsj/n5uwulphZZtSsOUKN0S2ymTmKEvaFv8YbGYOsBrami5XLrLODSxpjxqxENG+AjvZhJI+laQWGrDUXgH5sSRlhuvEmNePx6XTwFm98SvtZPkZ+1AAAgAElEQVTWA/pAaMEO4dIZxNQzmYXoSseyWYuxmXQYYW/oGdqufQ05+kU79RJgZCm33upV3TFpSAkzx7Kn6wsIW2Yn1GghKvDrhRGeKeKVC9WQz2L+NWR0TYLvmROaZJbhDmF4L3LzLtxSTYh2iybd7Sur24kR9oaeoCPal9YbRsLsqzDQ7nJd5MCmapz6uvZx6dWEZmXi1j57JPs9BIXXyH418bv1tlrn5nmZ7OGhkM9LZb030uyx2pyHLuQxhOYcojHm26lJN9yJrGBsptXGCHtD19Mp7Stuo/aUUD33zZZW8GrLlQJ23leZGFXCLX4MypZemZgUc6/WFVzaa2nMKuLyOSgXkMO3gpSIiycbbjtrsYynmUC1p/4Jr1xQpp/hcRh+S9w0VqmHHL41/iydO3Hda4iLr7X0PJtSCFYoNlM3YIS9ofvpoPalnRBsQ8wTv1w5fCuM3a3tqKLXBuLCPSGUsfSWleeOt4xYmkFWok5Kbxl2H9YKryy28HqacdIEqrU4jQdqAVZSfBbP1a8jEHasI2yUZhWCfvIAMsLe0P10WPvyhaDvKy6W2leuN+Qo23loh9JwPV9gzhRr++pFpHTuVMJRSmRl1amsRImU/ohhvdOU8MqqGUc7jWoHdN1wslnKH61cfA25/V79ca2M2FpQCFbLA2ilMcLe0PVk0b5amfDzhsaQO94Fm3cpQeoupbpeNkTSKtmxuxtK0K3cJ2t+8qFVrpA6YsjSHmmasb2UHPsm2kEk2uynv4N9/ll1TvBZRrFs9Sx++OW6dQ7RokLQLdmsOjlxbYS9oSdIE2CtTN6GzvWp/N8bGlWxZVogblN3lf2+QXODHL61vnkpoMk2LLzSNOMEYe9uORDvIDwPFcgsnKZQzJ9Qnerm8WqyEwZvCJu4fDbvqk1kZ2QtmGM67QJqhL2hZ9AJsFYmb2PnRhkcrQYRa4WQ7/rAZtj5zvABWcwNKQGuqrQ7Tk+kPG9oFDmoOls5sj880vCxLLh8DjbsrLlllp6veOAEEopLD64kzI+I5uZjWjHHrLYr6Eq4gBphb+htWpm8refrnmK+aJTgvIDcfk/D5gZx8WTF1h0RjsHMTjPHYMjBw58vyC7AUjVjP8TC+EMqSqc/OklYoyDO/YP6f2DiWe57OFx3YdWiXWrKaLbTCioEWe+/KxZVrYALqBH2ht6mCVtt0KMlzdfdWiw3Hb43iWbNDeq858PnXTwNs68oL5yNN8PWA8htb1MCa6GkIlg2IMCSNGPPz4hVXe1r11bVBpESZl6JeRJ5W3P6Nk4S9G0wv2QV4F2zqGoFXECNsDf0NI0Kz9iE4tIMDG2LxICpv4K2lWF/s+aGZGE8plwegwIrqDU3IMC0tv7BUaphjtMQAnHlR/HtSQvIYrF3PDh9pOl8tj4NCfAuWVS1EnMORtgbElltO2ZWUoVgYJtWCAxuhfnXYPiWgI35xXj+tQDtGPY36/2hPS9jaOIsAkz3zOXGG4HT9SuXEn5BCTKNzT4w+mglcXmIRgR4Fy2q6rQLqBH2hireyG3IxXLbwtauZGcRFYK6+oulGa0QEPMn4fzRal3tpXLiCtp2DPvb3S7SW467PEa15ooAS7u2ts1mJ1SohMt1hL2UkNIG1UVmm8fVCGD+hDot8Ltt70gDArzbvHg66QJqhL1B5TCljNz1bqS0EsPWNiLQVnPSK0kgy9NHEoVA6CNLc3Fscdgfa5fZYyqUQXD00UBHUC1PWDUBH/0b8ojRP5PENnOvIdImsX2EQFw+m3pIsI2j7SAtu21CrlEBbhZVGfoCb2gsHL/dspPD1mYUaC27Q7b60SUJ5MEb4NKZ2uKpBCHgT0hq/exbGPZr22Xr7cgbmptUjYdJFmGNXoiqHVy8OVPxiEl4Jklt5pu26t5cdtNH5hDQLbwHjQrwbllU1UmMsO93tPZeTdjaRuyYTWq/bRsNpKxarXmyvK6Nse7ufACx7Q6VYHv8IbwLL4fq0NKwP8223uCkqjc0htyWr2+rF5aKa5/wTOSuw3hnjiR2YmL+BFgCInnRw5XJFkGzSp33o52Zr7zK9Xx31H7GCPt+x//Ig0TD1jZqx2zWHbJNLnCZVq1u2gVTCXUQXmodmh32S8+N29KDZJxUDWe40mST0tjq5fVb9R4x60eVxl9+Qf3TdGIiZQ6D2ZOI6aONPaOU96Od70FX+M93EUbY9znW4jTehZdqmpsmbG2jQ+mmtN82u8A1tWrVr4P06tYhcTWvpr2q2wc2Jwt6SBTUQWIhCqKmGyR4kmCo5JCtXkdFoIrjj6lEJ0nPXHph/3rpVQV9NBpnQ2GZA++H1ie/ifega/znuwgj7A3YU8/B3n2IM9+EijcOtGbHbFj77YALXMOrVv06BAVaxjokaZFhLdxVJg8ryWsmLqhjMYB0IQqCHYSwABcmn0NcfE2VGrTVJ1ERqNZMMfmZR1eYSc29B6Jwpk0+i9mJ6kgnFFu/Xe9Bl/jPdxNG2BuqWLOvIINabQJZJ88a6Sw65QJXrWsGs1S1DtvuqJzsQrnQeIydihbpLpQjWrgNeDVhVg074E+ohgV1bHTge93UQ9iIa5eSNWXtTcQFarXtNt4InCIWqsGyVLKS6EjDb4OkyedIohZpr6u+J1oT3KUMPv5Rush/vltoSdgXiqWbgM8B9wMLwH/N55zfSTj2a8CPA76B+LV8znl7K9c3rDydtIO22wUutlp2dkK5B6aUbU8+jZybgO0jiJOPY6UEQquaLgY267XITbv1k9/nv4W4Oq83LwUEdYikOD6eqwRssBMICjWt0PMq54jab90ootJ2wrsKl07pr41I7kySJp+DydU15pWqT/6Od6m5leFbkJt2N/SudZv/fDfQqmb/Z8ArwAeBncCThWLp9XzO+R+aYw8Cd+dzTlGzz9ACK7V4qVN20Fj923AP2rqO7IfyCxnc8MrASGpsnLrmGc9VAnLrAa2Hi9WAeckbGlMdQzRpeMWtkoRkJYmjmqCgBxCoFIhJbaeNYaM6CJbeSJ90rl4jZX+SeWXTrpbetX7xn89K08K+UCzZwCLwyXzOuQqcKhRLfwu8E/gfkWO3AaPARCPXEEIg2pAiLlymFfrb67g3/wSM7Ku6S3oXXlI2+AbI2iZiaFR5qgRNPUJtF5EIkd7QqIqpslSuCM+E+u+4T/n5VwRRM/Vvta6xc+u0hzc0qtwzQ+V74AYShJdfxJ5/Fbe8LXR/lF9UK3SFhb1Uxi1/N3F/sH2EZYO3DGK55hpbuQbzr+LNTYTa2935rmq5wnNVwLTL55DXbYbtd8fj/+x4J+L1v9W2nZDXQn+r97tUUqMXljX2/DpeQqFjXcRSKdTerTy/IPZSuRa9tE3ffDfLEJHSqQqZJU52BgrF0gBQAP4on3M+F9n3buCvgCLwNuAl4N/nc86rsYKAT33+G8PAXOH5Z3HdDAs6DAaDwYBt2+QP3Q+w5ZGPPjgf3FdXsy8USx8E/lqz6wv5nPORyjE28OfAVeC/aY4dBL4NfAw4BXwc+LtCsXRbZVSgZfL0SZavXUva3RRCWOzcM87kqZOZJiO7FW9oFDn+kNZeKs58E2v2lcxlNdImUU2cS2cR09+qau/aenmusn9HNHxv5Dbkrnc3XP9mRw2UX8Q+fzT1/qB+ezRyj83WHVponzrnuXs/AJt3p5YbbDvhLnHzlX/k7KYfR4oBdXCS5h5o5+i9+r+lXA5H6QRt+8We30IZ1o/WfrdpFNgo3SxD1g0M+MI+vi/D+V8BNmm2XwMoFEsbgMcBB/iJfM5Zih6YzzlfA77m/y4US58Afh2l5ReSLiyl7FhjSul13YNqBDnoqEQSsR0eLJaburcsbWJNPoM3e7yWs3XLrcjNt+BVJs+09bIHYNBBRiY75WIZKS2NzTq5/u6u99XCOXhu9bppdQ3abBsZx7pb9lcDwwURC9N4F14O2+wvnwVJavtFJ4zT6g5++9iRuYD6z1fbrlIiN96MmPk+TB1Fbtqb2u7V57x5HISE9SBZVxP2UWuBAC58H3GhUG1nsTBdzfQV/C235jK9I8HnJ71l2H24Vmd7AEbvwp093j6PrQbt+t0oQ9IsNXWFfT7nuMBl3b5CsbQVeBKYBP5lPudcSTju3wB2Puf4WYRtYACIdQyGjGhDAkiY/2FbJqLqfgCbdtUmCwOTZ+2MOBitg7vrMGy9raE47c1M+MYCw/nRHwN1adRjpOnJbZ1Q1ZUdqJu1OK2C2QXbSggYOYBXcSXN4qlSXYwlPLj4RPokq+dWBX1dEt5dufFmmAn7b1TXSmzNqWBpoZ3x5ORNB5Lrg1W2rXrj/HdUoOt/W+kUkhgCfr9QLL1cOf53UV48bQhe3Z/E/ZErw+jNe3F3PtDSC1v3A0hZsGLNFNsScTAeHbISibPJOO1Z0QaGGz2EdO6sBk8LtUdWj5FmFvkMOZoFVFbonKRnJS6fRd5we+L16nmqhDonX3uVMrKCNhxVM6uSoTqjY7D19khntB8vyVsqKQFKIDl5o4K731bZtuKNcxB4P0o7ny8Uq5rbX+Vzzi8ViqU/BsjnnH+Xzzl/USiW9gL/AGwBngF+Op9z2p31ra+wJ59WC3f2HK4JhRZeWG9oFG/zOPiCLam8Otp7qxEH9dEhNZE4oWK2auNCGa1QDlw3OIrpdJKMOuekCass1wuuMPa25sIJS4Zv1awRELW4+RnXLSQhLp9TC66CpHR+1uI07qUzMLw3UpBdDXTWsODus1W2TQv7fM55Ce2gsrr/30V+Pwo82uz1DHqEZSNbCEccRI4/BLp5gEh5WcwALfnLZ43EKSXMHGuvFqYLDBfFspHDt6qIkB1MklH3nDaMsHTaMKCySkXxXDjzjbqxb0KnJJlVmuj8xNQzyM0B82HwnGYEd5+tsjXhEnodrf2zsRc2llA6dkC8vI4uWEn6CKOLg2aOYZ890r7rgj4wXHSFKsDY3cqG3OEkGanntDjCcrfcFo61Y9lqVOeHMogy+2o1baBuNBDF3XU49LyCZpWmO79SQgA1/94bENz9tsrWCPsep/rCBj9aKdTkWsJL6w2NhdLBibSE0ikfQKcSPtQmGGteN61G4myEaGA4ObI/bNqCqpmgbqRIzb012mZJ57QywgpluAoi7MTxurgyGT43xTauvKbSJ9Pb2fk1K7j7aZWtEfZrADE7gRy9MxBvxEq0V6oPtZb4WW6/RwmDa8+HC/VcmP5Oe3ODZsTd+UBlMtaqavT21D8BnetgdFQDwy1O47rLiWGSUyNFdphmhFXM1h9EuhXNXrNvqZxpUtPdclulo64/md7Ozq9Zwb2S79Rq0n3rfQ0NkZixyP+wIseqDzXw2IWl0vUFkRIunsI+/+yKC3p9TJsDavsqIi6ejNvyGzSXdQprcRprppj9WSUFVZMeTD+vTFOeF94OasFTmm2cSke9+736yfQVsIc33BZ9hNHse5jUjEW6DyvpI4+lshMw3LoLZ1N0qYdEM+ayriXJ1n/6SM0mPztRM/VdPAHbR9TK4euGE23j6SMGCTOvGCG8ihhhvwq0I0plarLpJHtlkq+ybtWdCJuConVu5h6C5wBt89JohUbuI8lc5i6UGvJQWW2S7Nu+oPeP8TsxFfBrpOatJT2qUT7rZZkCdfzMMeyzT6zQHRp0GGG/wrRtxV5S8uoL34dLpxDWOpXCbvCGUDYg9ZHXbPZID64kCKiKRu0G09rFElG4uFfOw9UFdf25iZDQ0N53JHlFq14azdLws0gadex5n3J/7aEVmI3Yt2PeWklJVjKMGAyrhxH2K0hbV+wlfVgC2H1YuQUGTDty+7245edry/wD3ji2n1A6mmPUc1VMkmDnYNmRRBR2JZNRhS1vwb1yJ/aJLybfd8bkFcGcpv4qyXaR9izspBC6SSEq2rCgreG6t8F7JPPEpM5bS5NkJcuIwbB6GGG/krTRHq39sGYnan7NELbFBzx0AMS1izVh4Qur8oswelfoQ2Vwa9znOi1GihCwYTvuTQ8i3nhJlZ80V5DSBtbidHVEIVsdBelIexYJwl7b5is8v6ALI9HsKtbMLKlkLiESTGv95MrYaxhhv5K02R4d/bAYcuLxUIJYdi1wV0CAWnOVDmDuOESiRLrb9eFSUxECRu9Abrsd9+IZmPmBXjD6JGRn6mjckiafRXTUEYrEmLGMZtGHkbhNPfOssWCaEMJ+5i48V0WbrGNa6xdXxl7DCPsVpBP26OCHpV1FGES6msBdh5Cjd8ClJ5HjDyEvvBwSGGL+hPLFD7rS+YGx6mXqEbaKZbJpF6QknNZGXNTFZmmD1hyaJE56FnXuK9jmbkJKwI6QNE8DdTvDdswViZOPw6DR2HsVI+zbQCMaUyeHuTX3wMhqT1C25asLcH0kNYFlgQxrz1HvErf0fK1M6cL084i5Vysx7XfXF/qWrQT96SPVcoHkZfy+YIoSDQLWYBvGEpCXX1ArYFt4FitqtkjypvJJ6AzbNUqyFsuxnASG3sEI+xZpRmPq5DC3OgG756eUYA+GkL1ufbqwALXPn+Ct3I+YnUB6ruow3pxRAhuwf/hltahr+Fa4fptyy9t8o174WzZsuBkWp9TPhDZI9NWWNa25mTZPEnjMTqgVsC3Q7PPM5IoauY42rHW1wAQTUpeuXfBp14SzIR0j7FtAK0CcQ7gLpRX1QHC3HIBNe+DSqdp1r9ugWSxlw+UfKW8aX7tNEsz+32As94r5RQZcJsXsBOLqPFTCKrg7H9CPLEDZ8cXb0wV00mTu/OuIuVeb11K7TOBldUWNEpoz2HhTODBckgmpi6M7rqXkId3eaRlh3wpJoXh3H8Zd76zIS+ve+nM1V8gbbse9cifije8lL27ZuL0muBfOK7fJtKwC0dAKIZfJSFKPimeIPHVEJaYIzg9ANvtyknvjlnHk5r1w6UxzQruDAq/Rj7xRV9Qo1ZHETDFTYLhm54r8+xJLKzjh3KPJQ3qh0zLCvhWSbKgr9NK6W26L+Lwrt0d5ZUqfEJpAqF5hqXN93/xmiCb1CHiGUH5BdSY73pFwrl5Ap5oqLFvNETQhtOPlunDpdDN3HaKpj7wJV1QdnZwrCt2Xu1TL3NVOumy01Sy90mmZQGgtYC1OK6GmSzrsv7SdZNNufVyb9WP67TqzzqUz+vpDxbzQQELlkOZ+CJbeiAcP80kR0Pbk03D6CFz6kabOFly9VCu3AQ8Ye/JpNSE7/7oazQzfgtz3sBJsTZD4kdcL2lYvQUqGzsvd+QBy38PI3e/NfA9Zg4Rp74vAStp2oWuHLjEvNUSd4HDdghH2LWJPPg2nvt7US+sngGgmoqM3NFYRxhEbjJRw8Yf6+sho1EYJy1fgzP/SXCDgXulfw4+JUi0vpSOwLLh+pBJB0a1dz/+7UPGq0dy/ipx4GDbdqI/bc90m5dlz+gn1983ZxtpQlzc2w/mx59XkR15VEqrtEmnXOp1X051MVpJGHoPtFfaxdmig427l22k7PdJpGTNOG7DnjjXsb500/M8yNA9NgkpZM3VICVemsEvfxl13faw+DGwMJ5QQQk3wIYFTtQt4Xngk4OcePfV1xJsz1frJYMwcHQObsM99Q8XSCS4+qow+5FsfjiXx1gZ4i2LZCGsdcnBr4ytsUwS0V9mvs1HrnpeYnWh6HiC2IM6vW+S5a9+HTps/ksyTSWEkWkDXDnUzYHWRfbz6fGYnVA6GlVhv0SRG2LeJRl7aJM3MXbc+MY2bf54cHg+H2fWF/OzJUBAynY3WveWD+oQSI/vg0qnaNq3Jx1IeP2/O1FwVF6eVm+f2+2DLW+KNsnylcuo6JZBD120giXcUz0Wu3wnbbk9PjK4jYaJWbrwJbn6P1kad5rbZSFrCKFGXTb+z8Sr7EoVah71rtBO6VPzs23IFPXL0zpDA1AnxbrKPuzvuQwbDi8weQ1w+Z7xx+gH/462reSRpZsG4NpGXOFRmFCFgqRxz9wytrh0aU+aLKFLjfqnTpqWEbW9Dbj0Qi1LpnX8OObw3tspWzJ+o/NeNTxhHsexqcLa4IPNAyMoowKuGY9CVUU+7zRRTKGCjFgsJsX0q12rHoipvaKyyQG1XyLOpKvj86wXfhw5HBg3dl7ymoiX47dFGEnMyJAnxbprU3XYw/HxGDkC50JWCHloU9oViaRtQAhYCm38nn3M+pTn2fcB/Bm4C/gn4cD7ndJdRqw1k0jy07oVeunkhzVzin59GksZ8eQrWb4uUpRHMkclX13OrKQvlyH5lCfJPkR5Mf7d6v8KykWmCHpQQ33FPzS3Uj5cuPaUxlV9Qo5rt9+h9+Cv1ll44OqPODJI5ptDgKCxM19WkW1kkp+3ELbuWfzdIQKjpOpl2+3n7wejEtjvgsgqn4UXCabRyzbomO50Q76Y1A93S6WSkVc3+IPBiPudo1rbXKBRL24G/BH4K+BbwWZTgf7jF63cfGTSPZO1yv/4lrmfa8DzExdf0uyofo/SWtR+J+NE/KDPOYGVbPQ0clBDe8Q7k2N0BDTSo1UvE3Ku13/WW+UupOgpfiAsL8GrnjOxHLl8BdzlZ0IMK2WzVXum0EVammEIVG7U2D4BoLUtVdeVx4mjNqivUQjF6OmDHrgpjUVEk0kabzVyz7nsdF+Irme+gLn5guODvLpuUDdIOYf+9DMf9G+BoPuc8BVAolj4O/KhQLP1KPudcbLEO3UVGzUNrU1++on2JtcKoXlYqNLFggsHIAucJP56954HdgNN9Bg20yqUzNTNFFCGIOftrFnDVXQ8QjZ2T0babxUYtZifUArJg/ZqwFVdNNtEFZ1GksgFnWSHbMTu2L4yDo8ak0WYz10zMD5D+XndNGOULL8VCgnerCQfaI+zfUiiWTgBDwF8BH8/nnKuR4/YDE/6PfM4pF4qlBWAcKCQVLoSopERrH3557S7Xx14q45a/W7PneS6UX1QJMSLXtJfKNQ8HYbFu6hncxWnYuBsun8aeOw7C0pc5d0JNBi2VVQjaSNne0KgafguvlpRkaATOfB0h1oXOk+tHK23iKvtso7hX453bdZvwRvYrYbXxJoRlg7cMwo2PHKRbCReQprUHjq2eFxiFeB6UC9V2FkOjtXsPlCGGRlXnFmHd1DN4cxMwOIp4swxjI6F3RAyNAtnL0+HuuA9GDyKEHS8riufB8mXEiceUOSnhOVfr1sC9ZkUslcBdQlQ0eyGvqdHgUglauKY3NKraeakM0fd69jjiyo9S7xfi385K4r8X685/C3fueN3ns7J1S9aIhNT5MQcoFEsfBP5as+sLwBVgFvh9YAvwZeCJfM757UgZfwqU8jnn44Ftk8BD+ZzzbLTgT33+G8PAXOH5Z3HdlMUnBoPBYKhi2zb5Q/cDbHnkow/OB/dl0ey/AmzSbL+WzzlvBn7PF4ql3wM+Afx25NgFlOYfZD1wOe3Ck6dPsnytMU3T1xqqPW0EISx27hln8tRJZCOrQzuMNzSqEjpHtVvpwaWziPNHtfcTLcO/dyBenuciTj5eLce/phAeN1/6e85u+nGkCNggpafXVNLs+lls/kEuTcGGgGnp8jllvtj1YP1ypr6tbKZBzfDCS9hTzwEVTTo6wjp/tG6V5PpRbhob4dz0LGKh1ubNlgfgjdyG3PXuhAu6MH8attwS2yXOfBNr9pW65bdSt3q0qz2073jknQwem/YdrybdKkMA1g0M+MI+vq/eyfmc46IRyoViSRSKpU8Cf5TPOWcrm68HljTFTADvD5zrABuBk2nXllI21JhRG7WXMmEkpddVD0oOOkh7ML5DAMPjyE17U+8nZp+fPQYXXo7Eby9gLUxX7dDVa1ZMN1IM1IS9dOHyedh0UzzGjpWWlrCRm3ZVIDb/47cHYNNexNRR5OyJ8AKwKF4lUqRvM/XPH70Ld/a4Cg0w+QxeJPNWPT9xd+cDVe8T7y3hZC7NlFe91cUyUlrxztxzoVRQ+QGGxzVzPeVM72krdauHEvAjiIVwXRq9pvYdtwdg0AnFyW/kO15Nuk2GgJKZSTRts8/nHFkolg4BjxaKpV8GxoCPA3+iOfxvgN8tFEsPAk8DjwJfzeecVM2+EbppsUVT7miLpYrPe4LdWnM/IU+bWLq622HmlfTkHP4EWSjJuFfzjNl0o94Fs1HtPVZ2xdPk0mkYjmizvv/6ma/jImHktpoHTHClsGXBjndqPZ/k5nG84D1X3BK9rblqAnP/b7Bd6nmfQPNultpAbBdPI6aerZbdSmTK4L2uJFUHgsBisEQyOC9003e81mh1gvYXgM8Bk8BVlKD/A6h63Lwrn3MO53POZKFYeoian/0zwIdbvHaYLlls0aw7mrU4rfK1Du9NOah2P7F46LqFUVsPQPmFxOQccmR/WGh7FV//YGRMXfTMrIK+Kpxrgg0IrTKWm3Ynfvz2mSO4iIDHjx8awqq1RyyBh/LXl4EwDEBo4Y6sdBiy4vVRfUYp3ifteIfqeZG0FJlylcIGNFKHTG6TXfIdr0VaEvb5nPMj4AMJ+z4Z+f0N4LZWrpdKFyy2aFkrmXkFNu9JMV2o+1HXCfp8a4RyZbscvjU9I5Swah4uQuo7Db9dG9XofUE/cwz7zJHa9kB9QmkUZfjj94bGlA9/sPPRXiPgridEOISCc6gyEogs3NGs1NSPdNr7DtXTvrNq592gATdTh7odWhd8x2uVNRP1spUIem2jxVCnwrLTBX3lfpRtVyOUdfa6sbv14W+1iVciWi1Uomieggvfb850I2wYORCKTuibVNwtB9QkXHUSgfA1ssTJUSXC1Ldg+juazkpjJ49SeUZd8Q5lpRvC6rYQ9TMp1HJPPYMeY03Fxln1xRataiVJ509/pxqaQJEgdBfKKvRBUOAlaVtJ17p0RtnRgxrw5j0qvvzWA+mCM0nzTzQ/RY6P1rXeylsfYavUiIsl5NjdGTuIAEHT0eTTyLkJ2D6iEqtf/GFjZa0U3aABd6gOq/4dr1HWjGbvk+g3OTMAAA61SURBVDVBQ6eu3YpWknS+ff7ZUBli/oRGA/cQZ5+AK5OaguPaVuxaAOUXVdm6GCWb9qiQDsHjo1yZ1O8PmZ/qhy/261pLDlNnrYV0a0Ihek+xY2X4b+QZuVsOIG96j9q3456Wkpt0km7QgOvVoZWY86v5Ha9V1pRm3w3U00rqeepk0WqsxWnc0vNhW/f082rnhh3xSnmeVtvyryWGRmEL2OeP4g6O6pewb3tbzaXzzYu1oGXVYzzVWbw5E4vgWA37sDVXX+uOaIb25NO47jLsfKf+eOnBdE3A+Pckh2+FqJYvPZh8DrF8ReuNU83nyzL4QTwy2sJXI9l0N2jASXXohsljQxgj7FcQd9f7ap4lKR9Alkk6bdTDrbmE2DPhn1HBJJbKsGVf9dqpOWBHDiCOP4a010WOUYnWZfkF7B9+WS/86pll/IBwETc+cfEkcvs9cXPBGz9AvPGSNtYNi9O4m/bUcvRKCVfOY5e+rb10KJ9vdOqjjjfIagq21XC3rFeHViaPV6PT7BeMsG8zSR++u+tweJFQG7wnYh/6YinBDdPS28w9F3f2GNbls0AtNHBVO96WVxp96KLKj10szSCnjoZ93WM2d03AsYVSWAAHXDNZXoCRA9Wk5X7bJbns2VP/lNg23tCYCvoWnHtYP4Y3NKZvb10+32phyXbobvCK6TqadJ80o4HOsuZs9qtJYgaqLQdUQLCkeN1twlqchpljca+cirBytxzQLr7yl/G7ez9Qta9ai9OIC4W4/dv3Y9/9Xth5n35R0/Ct2vppBbDnwqkjygffj/Do1y2QV9VPFi5OP4E4/liqoAca9xS5dErvzVTPFt4NXjHdRhM5WTueV9dghH1bSfrwN+3Rmy6k3pbeCvaZr8PMD2KTZnJkP+x5n8bdMtABbd6tJiRv+WBNA45OwEXjzrfq7mmpfLJZhGZDk3YNChx77hhcOR++n6X5+h1LjySbXkmamjw2nWbHMWacNlAvQQiXTsXdFqWEmWMdGerbZ47glQvhlapvfThb+FXLhuG9yM27cEsvhOYG5HXDsOMd4eOTPGoacPeU3rJKZN5GN75mklzYJ76oRj+bdsEA2K9+oW7sk65KptFFNJxIvBtcSdc4Rti3SJYEIfbcMdz1Tvi4mWPYZ4/Uv0CThLIY3fJB/cRt2opYEbe/e0Nj2f3YNTZa7eSvZVcndusl72508q4ZbxV77hhi/lWVzCUj3eAV041UJ8oz2OJNp9l5jLBvAa2dcb0Dp4/E3PqaFQhJAi6r4EtONO7B/A/rxuKR2/J4F1QS5VRPHV35Ce6e7kIZ9hwOx7gZvVMFbUtoo1ZiDq2Et0o3eMV0I41MYJtOs7MYYd8KKTZoXfCxRgVComdPI4IvKeTAxVPYr39FmS12v0d/bsW/Xm49UL2GEtYl2H04WcOvY6JSCcj1qQytmWKsjYzHSw/ToGeO6TQ7h5mgbYUOTs6levY04rWQUMdqBMqtt9dPUOJfe/v9eENjatSS6CvvwdxraqI4iUbbzUze9S5mArtrMMK+BTq6ZL0Rz54UwZdWx2QTj8Y8Y9mw4x3IfQ8jN96cHJLAsmDL3tQQAw23mxEYPUs3hHUwKIwZp0XaaWcM2uGTA5Wdinv2RASf8ijZA5dOYc8dS66jLn47qBg367frtXfLVmGHZ4+F/eKDRCd3NTTSbmbyrrcxtvjuwAj7NtAOO2PMqyfBO0Xr2RMM5OXHdxECbrgd98qd2Ce+qK9jQvx2ce7vlV9+0FwUumEbcfkc8uIp5buvMwMl2GWbzawkZieQ7jIIEYkAaugFjC1+9THCvgtIss8neackBp8KxncB9XfDdtwtB9SioQhVjXnbHZWKuCpPre9qOTuB3DweD3pWGUnYi9O4JUffKURGG97QmAqQtmlXwx410Y5QWrYRHAZDgxhh3w2kTEDqvFMgQVPSxXcRQpl0NMIewvHbxcnHsQKJn6t+0qGgZ5XcsYHzvdmJxEiXUBHWfoTO4P1ljSZpPHEMhpYxwr4baNfqwUun4Ibb4zljL51KPMUbGkMMjgLLWItlXI3/fkigb9oFw7cgN+0OBSojIdJlLf1hgv2/Xm5Rk5PUYGgLRth3Ae2agLTnjuFeuTMS1ndKa8KBgHlEeHDxCdy3PITcsDPZzOKbYKDmjrn8JmL5SrINPi21YJYOzSyjNxjaghH2XUK7PBZq8V32VL1xdITMI743zoaxZHNJkoa9835kJXql1gafFMM+Y4dmPHEMhvZghH0X0S6PBXvuGMwdq6aF03Ye2oTjCSGYF6eThXad+PzxEAsqdr2YejazwDauewZD67Qk7AvF0uXIpuuB1/M5JxZFqlAsfQz4JHA1sHlfPudokqYaWqVuSIUsybwD5pKY0NaRYEtvh7A2rnsGQ2u0JOzzOWej//9CsTQGFIDfSDj8IPCb+Zzz2VauaahPFg+WqvB2DukLkV7MXFKNi5PkW59iSzfC2mBYXdoZLuEPgb/N55ykoCgHge+18XqGJDLGkrEnn4ZTX9eHIjj1dW3SDnvuGJSer53jJ/swtnSDoatpi82+UCzdD/wYcEvC/kHgrcAjhWLpy8AU8PF8zvm7tHKFEIgsCTcawC+v3eV2E2KpBO5SzINFLJVi971u/lW8C1thEIS8VllY9SL2/KuJyU7WTT2DNzcBg6NIuYwQ62CpjLVYzpYgpcvph3ekEUx7hOnm9hBJIcfJIOwLxdIHgb/W7PpCPud8pPL/R4DP5HPOfEIxDnAU+CzwFPATwJcKxdJd+ZyjdxcBdu4ex3UTAm61yM494x0pt2u4/GR82/YRYERz8BsA3Hzp79XPITIm71gO/E0qu3dZ8+9Ig5j2CNON7WHbyXNwWTT7rwCbNNuvARSKpW3Ag8D/mVRAPuecAYJhEL9WKJaeAg4DicJ+8vRJlq9dy1DF7AhhsXPPOJOnTtZNOdeNuDvug20Ha26IF17CnnpOe6w3NAqDozWtO4GVaJOsdekGev0daTemPcJ0c3usGxggf+h+/b56J+dzjgtEvW6CvA94Lp9zzicdUCiW3g48mM85nwpsvh5YSru2lLJjjSml13UPqh7e0Bhy9K6aecYegNG7cGePa23lYmEaKuEPNGnBY3SqTaKeQV7GmDirTS++I53EtEeYbmwPKZO/9HbY7O8BvlXnmIvAJwrF0gTwNeBngHuBD7fh+v1DD4YOMLFtDIbuoB0zDHtQE64hCsXSw4Vi6QcA+ZzzOvAh4HeBS8BvAT+Zzzmx8wwp9GISD5NlymDoClrW7PM55/0J2x8DHgv8/irw1Vav18/0ZOgAE9vGYOgKTLiEHqPXQgf0ZAdlMKxBjLDvQXptNWqvdVAGw1rECHvDitBrHZTBsNYwwt7QVegSoBgMhtYxwr5FjHBqH3UjdRoMhqYxwr4FjHBqH8Yf32DoLN0XyadHSBROQ2OrW7FexfjjGwwdxQj7ZjHCqb304oIxg6GHMMK+WYxwaivW4jSUX6i1qfHHNxjairHZN4lZLNR+jD++wdA5jLBvASOc2o/xxzcYOoMR9i1ihJPBYOgFjM3eYDAY+gAj7A0Gg6EPMMLeYDAY+gAj7A0Gg6EP6OoJ2nXrBtpephAC27ZZNzCQmq+xnzBtEsa0RxjTHmG6uT3SZKbotsoCfOrz37gZOLPa9TAYDIYeZdcjH33wbHBDt2r254BdqETlBoPBYMjOZpQMDdGVmr3BYDAY2ouZoDUYDIY+wAh7g8Fg6AOMsDcYDIY+wAh7g8Fg6AO61RunYxSKpW3AC8CP53POyco2C/g08POAB3wmn3MeXb1arjyVdikBC4HNv5PPOZ9apSqtCoVi6R7gj4F9wPeAD+dzzonVrdXqUSiWPgZ8Erga2Lwvn3MmV6lKq0KhWPq3wH/I55x3VX7fCvw3IA+8Bnw0n3P+eRWrWJe+EvaFYule4M9Qbp1Bfg24F/WBjwBPFoql7+RzzpMrXMXV5CDwYj7n3LnaFVktCsXSIPAV4GPAl4H/CPw5cN8qVmu1OQj8Zj7nfHa1K7IaFIolG/gNVIf33cCuLwGPA/8K+DngrwvF0t58znHjpXQHfWPGKRRLOeBvgP9bs/tngf+czzlvVLT9zwG/tJL16wIOojTZfubHgPl8zvmLfM65CjwKvK1QLB1Y5XqtJv3+Xvwe8IHKXwAq78M48J/yOedaPuf8GXAJ+InVqWI2+kbYA6eBt+Rzzl9p9u0HJgK/TwD99oEfBA4UiqUThWLpXKFY+k+FYum61a7UChN6Dypa2uv037sAVEc6bwUeKRRLpUKx9L1CsfT+1a7XCvPpfM75F6j3wGc/8Fo+5ywHtnW9zFhTZpxCsfRB4K81u76QzzkfSTl1A7AY+L0ArG9j1bqCtPYB5oG/B34f2IIyY/x25V+/EH0PYI2+CxlxgKPAZ4GnUJrrlwrF0l35nHNsVWu2QuRzzpRmc0++J2tK2KPsrZs026/VOW8BGAr8Xg9cbleluojE9snnnDcDv+cLxdLvAZ+gv4R99D2Atfsu1CWfc84ADwQ2fa1QLD0FHAb6Qtgn0JPvyZoS9pVhdzMNPkHN+4LK/yeSD+9NktqnUCyJQrH0SeCP8jnHD550PbC0kvXrAiaAX/R/VCbnbmENvgtZKBRLbwcejHhk9eN7EWUCuKVQLNmBCdl9KC+urmVNCfsW+EuUXfJplOb7K8Cvr26VVo58zpGFYukQ8GihWPplYAz4OPAnq1uzFecp4IZCsfQR4C9Q3jgn8jmnL4U9KhDhJwrF0gTwNeBnUF5rH17VWq0y+ZzzSqFYOgX8VmUE/LMo0+c/rma96tFPE7RpfBZ4GngZZaP843zO+erqVmnF+QVgGJgE/hn4/4A/WNUarTD5nLMIvB/4VeAN4D3AQ6taqVUkn3NeBz4E/C7K2+S3gJ9MsGP3Gz+N8t66gFIMP5DPOV094jFRLw0Gg6EPMJq9wWAw9AFG2BsMBkMfYIS9wWAw9AFG2BsMBkMfYIS9wWAw9AFG2BsMBkMfYIS9wWAw9AFG2BsMBkMf8P8D2MWBaQZ8ABgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "dark"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# create some fake data\n",
    "Gaussion_nums=3\n",
    "data_nums=Gaussion_nums*200\n",
    "test_mean=[[7,4],[-4,-5],[4,-2]]\n",
    "test_cov=[[[2,0],[0,1]],[[3,0],[0,2]],[[7,0],[0,4]]]\n",
    "X1=np.random.multivariate_normal(mean=test_mean[0],cov=test_cov[0],size=200)\n",
    "X2=np.random.multivariate_normal(mean=test_mean[1],cov=test_cov[1],size=200)\n",
    "X3=np.random.multivariate_normal(mean=test_mean[2],cov=test_cov[2],size=200)\n",
    "X=np.r_[X1,X2,X3]\n",
    "np.random.shuffle(X)\n",
    "print(X)\n",
    "plt.scatter(X[:,0],X[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-24T04:34:48.011675Z",
     "start_time": "2020-03-24T04:34:47.997746Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.4299562   1.87520843]\n",
      " [-0.44817686  0.98836333]\n",
      " [ 0.26589353  1.51867961]]\n",
      "[array([[1., 0.],\n",
      "       [0., 1.]]), array([[1., 0.],\n",
      "       [0., 1.]]), array([[1., 0.],\n",
      "       [0., 1.]])]\n",
      "[0.4220036  0.84150093 0.92263666]\n"
     ]
    }
   ],
   "source": [
    "# initialize the parameters for three Gaussians\n",
    "mu=np.random.randn(3,2)\n",
    "sigma=[np.eye(2,2) for i in range(3)]\n",
    "phi=np.random.rand(3)\n",
    "print(mu)\n",
    "print(sigma)\n",
    "print(phi)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-24T04:34:48.798854Z",
     "start_time": "2020-03-24T04:34:48.791321Z"
    }
   },
   "outputs": [],
   "source": [
    "# compute p(x_i|z_i=j)\n",
    "def condi_Gaussian(x,z):\n",
    "    m=mu[z]\n",
    "    sig=sigma[z]\n",
    "    sig_inv=np.linalg.inv(sig)\n",
    "    return np.exp(-0.5*(x-m)@sig_inv@(x-m))/np.linalg.det(sig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-24T04:34:52.728353Z",
     "start_time": "2020-03-24T04:34:49.328531Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "step 0, loss=-12683.83720\n",
      "step 1, loss=-3945.26756\n",
      "step 2, loss=-3090.45931\n",
      "step 3, loss=-2848.59902\n",
      "step 4, loss=-2769.13439\n",
      "step 5, loss=-2725.38915\n",
      "step 6, loss=-2682.83000\n",
      "step 7, loss=-2638.12836\n",
      "step 8, loss=-2593.73131\n",
      "step 9, loss=-2552.54702\n",
      "step 10, loss=-2516.40599\n",
      "step 11, loss=-2486.13103\n",
      "step 12, loss=-2462.52494\n",
      "step 13, loss=-2446.11136\n",
      "step 14, loss=-2435.91855\n",
      "step 15, loss=-2430.00561\n",
      "step 16, loss=-2426.67257\n",
      "step 17, loss=-2424.81086\n",
      "step 18, loss=-2423.77250\n",
      "step 19, loss=-2423.19303\n",
      "step 20, loss=-2422.86972\n",
      "step 21, loss=-2422.68978\n",
      "step 22, loss=-2422.59019\n",
      "step 23, loss=-2422.53561\n",
      "step 24, loss=-2422.50615\n",
      "step 25, loss=-2422.49061\n",
      "step 26, loss=-2422.48272\n",
      "step 27, loss=-2422.47894\n",
      "step 28, loss=-2422.47732\n",
      "step 29, loss=-2422.47679\n"
     ]
    }
   ],
   "source": [
    "# EM iteration\n",
    "W=np.zeros((data_nums,Gaussion_nums))\n",
    "loss=0\n",
    "old_loss=0\n",
    "for s in range(100):\n",
    "    # E-step\n",
    "    loss=0\n",
    "    for i in range(data_nums):\n",
    "        tmp=np.zeros(Gaussion_nums)\n",
    "        for j in range(Gaussion_nums):\n",
    "            tmp[j]=condi_Gaussian(X[i],j)*phi[j]\n",
    "        normalization=np.sum(tmp)\n",
    "        for j in range(Gaussion_nums):\n",
    "            W[i][j]=tmp[j]/normalization\n",
    "            loss+=(W[i][j]*np.log(tmp[j]/W[i][j]))\n",
    "    print(f\"step {s}, loss={loss:.5f}\")\n",
    "    if abs(old_loss-loss)<1e-3:\n",
    "        break\n",
    "    old_loss=loss\n",
    "    # M-step\n",
    "    for j in range(Gaussion_nums):\n",
    "        phi[j]=np.sum(W[:,j])/data_nums\n",
    "        sig_tmp=np.zeros((2,2))\n",
    "        mu_tmp=np.zeros((2,))\n",
    "        for i in range(data_nums):\n",
    "            mu_tmp+=W[i,j]*X[i]\n",
    "            sig_tmp+=W[i,j]*(X[i]-mu[j])[:,np.newaxis]@(X[i]-mu[j][np.newaxis,:])\n",
    "        mu[j]=mu_tmp/(phi[j]*data_nums)\n",
    "        sigma[j]=sig_tmp/(phi[j]*data_nums)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-24T04:34:52.751202Z",
     "start_time": "2020-03-24T04:34:52.733644Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[7, 4], [-4, -5], [4, -2]]\n",
      "[[ 6.89263494  3.87541387]\n",
      " [-3.79494432 -4.90925045]\n",
      " [ 4.27670931 -1.92711354]]\n",
      "[[[2, 0], [0, 1]], [[3, 0], [0, 2]], [[7, 0], [0, 4]]]\n",
      "[array([[2.12271919, 0.01537255],\n",
      "       [0.01537255, 1.50088058]]), array([[5.26666163, 0.4499918 ],\n",
      "       [0.4499918 , 1.69574572]]), array([[ 6.5063268 , -0.90276312],\n",
      "       [-0.90276312,  2.77246074]])]\n"
     ]
    }
   ],
   "source": [
    "print(test_mean)\n",
    "print(mu)\n",
    "print(test_cov)\n",
    "print(sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-24T04:34:56.562989Z",
     "start_time": "2020-03-24T04:34:56.541132Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "GaussianMixture(covariance_type='full', init_params='kmeans', max_iter=100,\n",
       "                means_init=None, n_components=3, n_init=1, precisions_init=None,\n",
       "                random_state=0, reg_covar=1e-06, tol=0.001, verbose=0,\n",
       "                verbose_interval=10, warm_start=False, weights_init=None)"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.mixture import GaussianMixture\n",
    "gmm=GaussianMixture(n_components=3,covariance_type='full', random_state=0)\n",
    "gmm.fit(X)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-03-24T04:34:57.496286Z",
     "start_time": "2020-03-24T04:34:57.489150Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 6.92629599  3.96447159]\n",
      " [-3.99933948 -4.97301208]\n",
      " [ 4.12592937 -1.90374852]]\n",
      "[[[ 2.06047734 -0.0538924 ]\n",
      "  [-0.0538924   1.27515566]]\n",
      "\n",
      " [[ 4.47676309  0.24102403]\n",
      "  [ 0.24102403  1.59594964]]\n",
      "\n",
      " [[ 7.0414659  -0.37232227]\n",
      "  [-0.37232227  3.16852078]]]\n"
     ]
    }
   ],
   "source": [
    "print(gmm.means_)\n",
    "print(gmm.covariances_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
